#!/usr/bin/Rscript

#plotAdcircOutput.R:
#
# Last modified 13 August 2008
# Author: Cameron Bracken
#         Student Researcher
#         NOAA National Weather Service
#         300 Startare Drive
#         Eureka, CA 95501
#         cameron.bracken@noaa.gov
#
#Overview:
#   This script reads in various files related to adcirc and generates a multipage pdf 
#   image which is then converted to a flash movie with the program pdf2swf which is part 
#   of swftools (http://www.swftools.org/).
#
#Inputs:   
#   This script calls the fortran program ppCurrents which post processes 
#   adcirc output and then read the postprocessed output files:
#      fort.64.chop  - contains the velocity data 
#      fort.64.times - contains the adcirc timestamps and timesteps
#      fort.14.chop - contains the nodes corresponding to the data
#   
#   If bathyContours = T then three additional files are read:
#      bath.cg3.sparse.grid - contains the table of bathymetry data 
#      lat.cg3.sparse.grid  -    "      "  latitides of the bathymetry data
#      lon.cg3.sparse.grid  -    "      "  longitudes "  "       "       "
#   these files are used to plot bathymetry contours 
#   
#   Also, one last file is read:
#      sfcoastline - containd the san francisco coastline information 
#   
#Outputs:
#   Using all of this information, plots of the adcirc output are created at every timestep 
#   available.  Two files are created: outfile.pdf and outfile.swf where outfile is specified 
#   in the plotting options section below
#
#NOTE: The parameter thin was introduces because pdf2swf cannnot render all of the arrows if 
#      if the output is not thinned  by arbitrarily taking every thin'th data point 
#   


######### plotting Options ##########
bathyContours = F                   #if T then interpolate and plot bathymetry contours
thin = 10                           #thin the number of arrows by a factor 
outdir = 'plots/adc'          #name of the output file without any extension
#####################################

con = file('/dev/null','w')
 sink(con,type='message')
     library(matlab)
     library(sfsmisc)
     library(methods)
 sink(NULL)
 sink(con,type='output')
     library(spam)
     library(fields)
 sink(NULL)

### Input files ###################
coastfile = 'chopped/sfcoastline'
chopfile = 'chopped/fort.64.chop'
timefile = 'chopped/fort.64.times'
nodefile = 'chopped/fort.14.chop'
bathfile = 'chopped/bath.cg3.sparse.grid'
bathlat = 'chopped/lat.cg3.sparse.grid'
bathlon = 'chopped/lon.cg3.sparse.grid'
####################################

if(bathyContours){
    bath = as.matrix(read.table(bathfile))
    xb = scan(bathlat)
    yb = scan(bathlon)
}

cat('Reading Data...\n')
parms = scan(chopfile, nlines = 1, skip = 1)
nsets = parms[1]
nnodes = parms[2]

nbound = scan(coastfile,skip=1,n=1,quiet=T)
totalskip = 2  #initially this is the number of header lines on the coastfile
nbpts = 0

    #read in all the coastlines (though we will only use the first one for now)
for(b in 1:nbound){
    nbpts[b] = scan(coastfile,skip=totalskip,n=1,quiet=T)
    totalskip = totalskip + 1
    thisbound = matrix(scan(coastfile,skip=totalskip,nlines=nbpts[b],quiet=T),ncol=2,byrow=T)
    totalskip = totalskip + nbpts[b]
    if(b==1){
        allbounds = list(thisbound)
    }else{
        allbounds = c(allbounds,list(thisbound))
    }
}

xy = matrix(scan(nodefile, nlines = nnodes),ncol=3,byrow=T)
uv = matrix(scan(chopfile, skip = 2),ncol=3,byrow=T)[,2:3]
times = matrix(scan(timefile),ncol=2,byrow=T)[,1]
times = as.POSIXlt(times,tz='UTC',origin='2004-01-01')

    #this little dude is the thinning vector (xy-skip-sequence)
xyskseq = seq(1,nnodes,by=thin)

cat('Creating plots...\n') 
pb <- txtProgressBar(min=1,max=nsets,style=3)
   
     
    outfile = paste(outdir,times[1]$year+1900,times[1]$mon+1,times[1]$mday,sep='')
    pdf(file=paste(outfile,'.pdf',sep=''))

        for(i in 1:nsets){
            ind1 = (i -1) * nnodes + 1
            ind2 = i * nnodes
            uvskseq = seq(ind1,ind2,by=thin)

                plot(xy[xyskseq,],type='n',xlab='',ylab='',bg='transparent',las=1)

                xpos = .5*abs(max(xy[xyskseq,1])-min(xy[xyskseq,1]))+min(xy[xyskseq,1])
                lab = format(times[i],'%Y-%m-%d %H:%M')
                text(x=xpos,y=max(xy[xyskseq,2]),lab)

                arrow.plot(xy[xyskseq,1],xy[xyskseq,2],uv[uvskseq,1],uv[uvskseq,2],xpd=F,
                                        arrowfun=p.arrows,size=.2,fill='steelblue')

                lines(allbounds[[1]],lwd=2)
                lines(rbind(allbounds[[14]],allbounds[[14]][1,]),lwd=2)
                lines(rbind(allbounds[[15]],allbounds[[15]][1,]),lwd=2)

                if(bathyContours)contour(yb,xb,bath,add=T,levels=c(-20,-30,-40))

                if(times[i]$hour==23 & times[i]$min==45 & i<nsets){
                    dev.off()
                    system(paste('/opt/local/bin/pdf2swf -s framerate=1 -l -B plots/./viewer.swf -o ',outfile,'.swf ',outfile,'.pdf',sep=''))
                    outfile = paste(outdir,times[i+1]$year+1900,times[i+1]$mon+1,times[i+1]$mday,sep='')
                    pdf(file=paste(outfile,'.pdf',sep=''))

                }
                
            setTxtProgressBar(pb, i)
        }
    close(pb)
dev.off()

    #Finally create the flash file
system(paste('/opt/local/bin/pdf2swf -l -B plots/./viewer.swf -s framerate=.1 -o ',outfile,'.swf ',outfile,'.pdf',sep=''))
    
